Hybrid Reynolds- Averaged / Large Eddy Simulation of a 
Cavity Flameholder; Assessment of Modeling Sensitivities 


Steady-state and scale-resolving simulations have been performed for flow in and around a model scramjet 
combustor flameholder. The cases simulated corresponded to those used to examine this flowfleld experimen- 
tally using particle image velocimetry. A variety of turbulence models were used for the steady-state Reynolds- 


eddy simulation everywhere except in the inner portion (log layer and below) of the boundary layer. Hence, 
this formulation can be regarded as a wall-modeled large eddy simulation. This effort was undertaken to 
formally assess the performance of the hybrid Reynolds-averaged / large eddy simulation modeling approach 
in a flowfleld of interest to the scramjet research community. The numerical errors were quantified for both 
the steady-state and scale-resolving simulations prior to making any claims of predictive accuracy relative 
to the measurements. The steady-state Reynolds-averaged results showed a high degree of variability when 
comparing the predictions obtained from each turbulence model, with the non-linear eddy viscosity model (an 
explicit algebraic stress model) providing the most accurate prediction of the measured values. The hybrid 
Reynolds-averaged / large eddy simulation results were carefully scrutinized to ensure that even the coarsest 
grid had an acceptable level of resolution for large eddy simulation, and that the time-averaged statistics were 
acceptably accurate. The autocorrelation and its Fourier transform were the primary tools used for this as- 
sessment. The statistics extracted from the hybrid simulation strategy proved to be more accurate than the 
Reynolds-averaged results obtained using the linear eddy viscosity models. However, there was no predictive 
improvement noted over the results obtained from the explicit Reynolds stress model. Fortunately, the numer- 
ical error assessment at most of the axial stations used to compare with measurements clearly indicated that 
the scale-resolving simulations were improving {i.e. approaching the measured values) as the grid was refined. 
Hence, unlike a Reynolds-averaged simulation, the hybrid approach provides a mechanism to the end-user for 
reducing model-form errors. 
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averaged simulations which included both linear and non-linear eddy viscosity models. The scale-resolving 
simulations used a hybrid Reynolds-averaged / large eddy simulation strategy that is designed to be a large 
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speed of sound 

user-specified constants for the Larsson sensor 

constant for the turbulent viscosity 

distance to the nearest solid surface 

non-dimensional “law of the wall” coordinate 

hybrid RAS/LES blending function 

factor of safety used to assess grid convergence 

generic functional used to assess grid convergence 

turbulence kinetic energy 

mass fiow rate 

pressure 

numerical scheme order of accuracy used to assess grid convergence 
turbulent Prandtl number 
autocorrelation function 

grid refinement factor used to assess grid convergence 
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5 Ct turbulent Schmidt number 

T temperature 

u, V, w Cartesian velocity components 

X, y, z Cartesian coordinates 

Vol grid cell volume 

y velocity vector 

a hybrid RAS/LES blending function constant 

6 boundary layer thickness 

A SGS filter width or grid spacing 

'^x,y,z or z direction grid spacings averaged over the computational domain 

6 small non-zero fioating point value 

T] length scale ratio controlling the hybrid RAS/LES blending function or Kolmogorov length scale 

6 argument associated with the Lars son sensor 

K von Karman constant or MUSCL parameter 

A turbulence length scale (Taylor microscale) 

Id molecular viscosity 

jdsgs sub-grid viscosity 

p density 

T time difference relative to the initial time for the autocorrelation function 

if/ sensor function used to blend non-dissipative and dissipative operators 

6g) specific turbulence dissipation rate 

V gradient operator 

Subscripts 

/, j computational coordinate indices 

oo freestream or reference value 


Introduction 

Reynolds-averaged Computational Eluid Dynamics (CED) codes are the standard high-fidelity numerical tools 
utilized in the aerospace industry. These tools have revolutionized the research and development practices, which as 
early as 15-20 years ago relied almost exclusively on extensive wind tunnel testing. Unfortunately, Reynolds -Averaged 
Simulations (RAS), which attempt to model all of the scales present in turbulent flows, have proven to be deficient in 
many challenging areas of interest to the aerospace community. Some examples include: 

• high lift devices (massive flow separation) 

• combusting flows (particularly lean or rich flames near extinction) 

• unsteady flows (rotorcraft, aeroacoustics, etc.) 

• shock / boundary layer and shock / jet interactions 

In general, the limitations associated with the turbulence closure models are the pacing items preventing the use of 
RAS in these (and other complex) settings as a true predictive tool.^ 

Large Eddy Simulation (EES) methods have the potential to reduce the modeling uncertainty inherent to RAS 
approaches, since the intent of EES is to resolve the large scale turbulent structures while modeling only the small- 
est scales. However, the computational expense of wall-resolved EES at Reynolds numbers relevant to engineering 
problems of interest is well beyond what can be deemed as practical. Hybrid RAS/LES approaches'"^ have emerged 
to address this issue, and have provided a rational path forward towards extending EES into practical settings. These 
methodologies allow EES content to be resolved in areas that require a rigorous modeling approach, while maintaining 
a more cost effective RAS approach for benign regions of the flow {e.g. attached boundary layers). 

The computational expense required for hybrid RAS/LES, while less than that of a full EES, is still formidable 
when compared with steady- state RAS. Moreover, the numerical algorithms required to resolve the turbulence scales 
of interest must have low numerical dissipation with minimal dispersive errors (particularly for high-speed flows where 
shock waves may be present). These observations demand an efficient high-order, low-dissipation numerical frame- 
work; a feature not typically required by pure RAS solvers. As a result, the extension of scale-resolving simulation 
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approaches (such as hybrid RAS/LES) to engineering problems of interest will only become practical when substan- 
tial advancements to both the numerical and physical models have been realized. Towards this end, researchers in the 
Hypersonic Airbreathing Propulsion Branch (HAPB) and Computational AeroSciences Branch (CASB) of the NASA 
Langley Research Center have been extending the capabilities of the VULCAN-CFD flow solver^’ ^ for scale-resolving 
simulations. Advancements made to the numerical framework^ and the hybrid RAS/LES framework^ have been doc- 
umented for both low and high-speed benchmark flows. The present effort aims to apply this framework to a flow of 
engineering interest. In particular, this paper describes the results of hybrid RAS/LES performed for a model scramjet 
combustor flameholder. The flowpath involves a supersonic internal flow passing over a recessed cavity that is inter- 
nally fueled by an array of ethylene-fueled injection ports. Experimental data available for this configuration include 
two components of velocity (via particle image velocimetry)^^ and reacting scalar data (via UV Raman Scattering). 
The results described in this paper focus (to a large extent) on the sensitivity of the predictions to numerical modeling 
choices utilized by the hybrid RAS/LES approach. Comparisons are also made with measured velocities within the 
cavity flameholder and with steady-state RAS results. 


Geometry Description and Flow Conditions 

A schematic of the facility flowpath considered in this study is shown in Fig. 1. An asymmetric facility nozzle 
provides a continuous flow of Mach 2 (nominal) air to the constant area isolator section of the flowpath. The 7 inch 
long isolator starts at the v = 0 inch station and has a 2 inch height (y-direction) and a 6 inch width (z-direction). 
The 2.5° divergent portion of the lower wall initiates at the exit of the isolator (x = 1 inch station), and the cavity 
flameholder is located 3 inches downstream of this location. The cavity spans the entire width of the flowpath. The 
depth of the cavity is 0.65 inches and the length of the cavity floor is 1.815 inches. The aft wall of the cavity has a 
shallow angle (22.5° relative to the cavity floor) and houses an array of 1 1 evenly distributed fuel injection ports. The 
diameter of each injection port is 0.0078 inches, and the centerline of each port intersects the angled aft cavity surface 
at the v,y coordinate values of 12.0276 and 0.7703 inches, respectively. The cavity closeout location corresponds to 
the intersection of the cavity aft wall with the 2.5° facility lower wall surface. 



Table 1 : Facility Test Conditions 


Nominal Conditions 

Case 1 

Case 2 

Case 3 

Air Mach Number 

2.0^ 

2.0^ 

2.0^ 

Air To [K] 

589.0 

589.0 

589.0 

Air Po [kPa] 

483.0 

483.0 

483.0 

Fuel Flow Rate [SLPM] 

0.0 

56.0 

99.0 

Fuel To [K] 

N/A 

310.0* 

310.0* 


^ nozzle design Mach number 
^ estimated bottle temperature 


The nominal flow conditions considered in the experiments conducted by Tuttle et al.^^ are given in Table 1. 
The fuel flow rates, reported in Standard Liters Per Minute (SLPM), are based on a reference state of 273 K and 
1 atmosphere. The facility nozzle geometry was included in the simulated flowpath, so the Mach 2.0 value listed in 
Table 1 is a nominal value. All of the simulations performed in this work ignored any facility side- wall influences. 
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This allowed the exit conditions from an a priori two-dimensional RAS of the facility nozzle (and 4.2 inches of the 
constant area isolator) to be used as the inflow condition for the three-dimensional simulations performed for the 
region of interest further downstream. The Mach number and pressure profiles extracted at this interface (see Fig. 1) 
are shown in Fig. 2. The simulation of the facility nozzle shows a weak shock structure propagating through the 
isolator as evidenced by the pressure variability along the profile. The boundary layer thickness at this station is 
0.25 inches (1/4 of the duct half height) which is nearly 40% of the cavity depth. 


Pressure [kPa] 


60.00 61.25 62.50 63.75 65.00 66.25 67.50 68.75 70.00 



Mach Number 


Figure 2: Facility Mach number and pressure profiles at the v = 4.2 inch station used to prescribe the 3-D RAS and 
hybrid RAS/LES inflow conditions 


Computational Grid Description 

As mentioned previously, the simulations performed for this effort neglected any facility side- wall influences. This 
allowed the consideration of only a fraction of the facility width, providing for a more efficient means of examining 
numerical modeling sensitivities while retaining the salient flow features of interest. The smallest spanwise width 
allowed by the symmetry present in the geometry is a 0.25 inch section that covers the region from the centerplane of 
an injection port to the gap centerplane between two injector ports. This domain is appropriate for steady-state RAS, 
with symmetry conditions imposed at each bounding plane. However, symmetry conditions are not appropriate for the 
unsteady hybrid RAS/LES, since the flow symmetry only applies to the time-averaged flowfield statistics. Hence, these 
simulations (at a minimum) require a domain that extends from the centerplane of one injector port to the centerplane 
of an adjacent port (or equivalently one that is bounded by the gap centerplane on either side of an injection port) 
to allow the specification of periodic conditions at each bounding plane. However, as pointed out previously, the 
approach boundary layer has a thickness of 0.25 inches which will eventually transition to an even thicker free shear 
layer over the cavity flameholder. In order to properly accommodate the formation of resolved eddy structures, the 
largest of which are comparable to the width of the viscous shear layer, a domain larger than 0.5 inches should be 
considered. More specifically stated, the width should be large enough to ensure that the largest eddy structures 
become decorrelated when separated by a distance larger than the half- width of the periodic domain. This prompted 
the consideration of a 1-inch domain width that encompasses two full injection ports. This domain matched that used 
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in previous DES research efforts^^’^^ for this configuration. 

A sequence of three progressively refined structured grids was created for the simulations. The refinement factor 
used was a factor of 1.5 in each coordinate direction. The grid attributes closely resembled those utilized for the DES 
studies performed in Refs. 12 and 13. However, the computational grids generated for this study were fully structured, 
while the grids in the aforementioned references made use of unstructured hexahedral cells in some regions of the flow. 
Each grid is comprised of nearly isotropic cells (except in the inner part of the boundary layers where RAS is utilized), 
and the grid is clustered to all surfaces such that is approximately unity. The first attribute is driven by the desire 
to resolve eddy structures in the outer portion of the boundary layers and in free shear regions of the fiow. The second 
is a resolution requirement for accurate RAS modeling of boundary layer fiows without resorting to the use of wall 
functions. The average grid spacing and wall distance (normalized by the incoming boundary layer thickness), average 
and maximum and the total number of cells associated with each grid generated for the RAS domain are listed in 
Table 2. The hybrid RAS/LES grids were obtained by duplicating and mirroring the RAS grids which results in a cell 
count that is 4 times larger than those listed in Table 2. Several images of the grid are shown in Eig. 3. Topological 
nesting was used in several areas to help maintain a nearly constant grid spacing as the cross-sectional area of the 
domain varied, and a non-trivial topology was implemented in the vicinity of the fuel injectors to provide additional 
grid resolution for the smaller fiow structures present when fuel is injected into the cavity. Also shown in this figure 
is the location of the recycling plane used to generate resolved turbulent content for the infiow of the scale-resolving 
simulations. 


Table 2: Grid Attributes (RAS domain) 


Attribute 

Coarse 

Medium 

Line 

Al S 

0.063 

0.042 

0.028 

di / s 

0.0004 

0.0004 

0.0004 

(ave, max) 

0.50, 0.85 

0.50, 0.92 

0.50, 0.99 

Total Cell Count 

2,540,336 

7,867,136 

24,585,120 


Uncertainty Quantification 


Estimating model-form uncertainties is a challenging task requiring extensive validation efforts at conditions that 
are representative of the problem of interest. Eor many turbulent fiows of engineering interest, the physical models 
chosen for the CED simulation are often the dominant source of uncertainty. This is particularly true for turbulence 
closure models in a RAS framework. This uncertainty can (presumably) be reduced by resorting to a scale-resolving 
simulation approach like hybrid RAS/LES. Prior to assessing the model-form uncertainty, the numerical errors {i.e. 
grid convergence) should be quantified to ensure that the uncertainties associated with the numerical treatment are 
sufficiently small. This error source can be quantified using the Grid Convergence Index (GCI).^^ The GCI is a grid 
convergence estimator derived from the generalized Richardson extrapolation formula and can be written as follows: 


GCI = Fs 


I/1-/2I 

rP - 1 


( 1 ) 


where / is some fiow parameter of interest evaluated at two different grid resolutions (/i and / 2 ), r is the grid refinement 
ratio, and p is order of accuracy of the numerical scheme. Einally, Fg is a safety factor with recommended values taken 
to be either 3 (if the observed order of accuracy is assumed to be the theoretical value) or 1.25 (if the observed order 
of accuracy has been rigorously determined). It should be emphasized that the GCI, while based on Richardson 
Extrapolation, is not meant to be a “best estimate” of the numerical error. Instead, the intent is to provide a reasonable 
bound on the discretization error. Note that the impact of aleatoric (or random) uncertainties should also be considered 
prior to assessing the predictive accuracy of the relevant physical sub-models. However, the information required to 
formally perform this assessment was not provided in the documents that described the experiment. 

The quantification of uncertainty is an even more daunting task for scale-resolving simulations due in part to the 
massive computational resources required to conduct each simulation. The level of statistical convergence must now 
be considered as a source of uncertainty, since the outputs of interest from an engineering perspective tend to involve 
time averages of the turbulent fiow. The sub-grid models associated with EES, which explicitly vary with grid spacing, 
also blur the distinction between epistemic model-form errors and numerical resolution errors. These complications 
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Figure 3: Grid images: overall view of the coarse grid (top 2 images), cavity fuel injection region (lower image) 


6 of 26 


American Institute of Aeronautics and Astronautics 



make any rigorous attempt at quantifying the numerical errors with LES or hybrid RAS/LES an arduous task that 
is seldom addressed in practice. In an attempt to simplify this process, the following viewpoint was adopted in this 
effort. The LES equations can be be written as the “under-resolved” Direct Numerical Simulation (DNS) equations 
plus (potentially) extra terms. The extra terms represent the modeled sub-grid contributions, which under normal 
circumstances, are explicitly dependent on the grid size. Eor instance, a typical molecular diffusion term in the DNS 
equation set 

duj 

Diffusion = [i - — (2) 


is replaced with a term like 


Diffusion = 


/ X OUi 


( 3 ) 


in the modeled LES equations. The Sub-Grid Scale (SGS) viscosity, [isgs, varies like (Av)^, so as the grid is refined 
towards that required by DNS resolution, the extra SGS terms approach zero like (Av)^. Based on this observation, 
one can choose to lump the SGS modeling error with the numerical error, which can in principle be estimated using 
measures like the GCI. In other words, the SGS terms can be decomposed as 


SGS = SGSexact 


-|- SGSe 


( 4 ) 


where the SGSexact is the SGS effect required to recover the DNS statistics {i.e. the ideal SGS model), and the SGSgrror 
can be considered as part of the numerical error. Note that an Implicit LES (ILLS) utilizes precisely the same equation 
set as DNS. Hence, the SGS model error for an ILLS (given the perspective described here) would be the leading error 
term of the numerical scheme. 

While this viewpoint on uncertainty quantification for scale-resolving simulations is convenient, the practitioner 
must appreciate the limitations of this approach. The GCI is an error estimator derived from the Richardson extrap- 
olation formula. This formula assumes that the dominant uncertainty in the solution is due to the smallest power of 
the grid spacing associated with the error of the numerical approach. This power is the minimum of 2 (assuming an 
explicit SGS model is used) or the leading error term of the truncated Taylor series associated with the numerical 
scheme. When this condition is satisfied, the numerical solution is said to be in its asymptotic convergence range, 
and the GCI can reliably be used to provide uncertainty bounds on the numerical error. However, if the grid is coarse 
enough such that higher-order terms are dominant, then the GCI estimate in general can not be relied upon to truly 
provide error bounds. However, from a practical viewpoint this scenario is often easy to detect. Eor instance, when the 
differences in solution estimates between grid levels are not monotonically reduced as the grid is refined, or when the 
difference in solution estimates are changing sign as the grid is refined. Hence, the informed user can often recognize 
this occurrence and tag the GCI value as unreliable, or add an additional safety factor if desired. At a minimum this 
approach at least provides a consistent means of documenting grid-related errors using a traceable framework. 


Reynolds-Averaged Results 

Reynolds-averaged simulations were performed for a variety of turbulence models to establish a baseline level 
of predictive accuracy for this flowfield based on the current state-of-the-art practices used for engineering purposes. 
These simulations were advanced in pseudo-time using an incomplete LU factorization scheme (with planar relax- 
ation)^^ with a Courant-Eriedrichs-Lewy (CEL) number of 100. The inviscid fluxes were evaluated using the Low- 
Dissipation Elux Split Scheme of Edwards with cell interface variable reconstruction achieved via the /c=l/3 Mono- 
tone Upstream-centered Scheme for Conservation Laws. The van Leer flux limiter was utilized to avoid spurious 
oscillations during this reconstruction process. The viscous fluxes were evaluated using 2"^-order accurate central 
differences with the constituent viscosities and conductivities computed from the polynomial fits of McBride. The 
turbulence models considered were the Menter-BSL k-oj model,^^ Menter-SST k-oj model,^^ and the Gatski-Rumsey 
Explicit Algebraic Stress (Gatski-EAS) model.^^ The turbulent Prandtl {Pvt) and Schmidt {S Ct) numbers, which con- 
trol the turbulent transport of energy and mass, were chosen as 0.9 and 0.5, respectively. Only simulations of Case 1 
(no fuel injection) are discussed in this effort. Reacting and non-reacting simulations of the cases with fuel injection 
will be the focal point of follow-on investigations of this flowfield. 

Solution convergence of the steady-state RAS was monitored by assessing the L 2 norm of the equation set residual 
error, the mass flow error, and the surface friction force time histories. A sample convergence history is given in 
Eig. 4, which shows the normalized mass flow rate error and the L 2 norm of the steady-state residual error for the 
coarse grid Gatski-EAS simulation. At a minimum, the following iterative convergence statements were satisfied for 
each simulation: 
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• The L 2 norm of the steady- state residual error was reduced by at least 4 orders of magnitude. 

• The surface friction force time history remained unchanged to 5 digits over the final 2500 iteration cycles. 

• The relative mass fiow rate error, |mout - ^inl / was less than 2.5 x 10“^ 
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Figure 4: Typical steady- state iterative convergence history 


Figure 5 shows the structure of the fiow (via streamlines) inside the cavity fiameholder as predicted by each of 
the turbulence models considered. Also shown are the streamlines extracted from time-averaged Particle Image Ve- 
locimetry (PIV) measurements taken by Tuttle et al}^ The predicted fiow structure is qualitatively similar to the 
measurements with a dominant clockwise rotating recirculation zone and a smaller counter-clockwise rotating recir- 
culation zone adjacent to the front wall of the cavity. From a fiameholding perspective, the primary recirculation zone 
provides the mass exchange between the core fiow and the hot cavity combustion products. The smaller recirculation 
zone simply enhances the cavity fiow residence time. The most notable difference between the model predictions is 
the size and shape of the secondary recirculation zone associated with the Gatski-EAS model, which is larger and more 
elongated than those predicted by the Menter models. This difference is likely due to Reynolds stress anisotropies that 
FAS models (or more generally, non-linear eddy viscosity models) are capable of capturing. The impact of captur- 
ing this feature will become evident in the subsequent discussions of profile properties extracted at various locations 
within the cavity. Figure 6 shows the specific locations where fiowfield profiles have been extracted. 



Figure 5: Streamlines extracted from measurements (bottom-right) and streamlines extracted from simulations: 
Menter-BSL (top left), Menter-SST (top-right), and Gatski-EAS (bottom-left) 
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Figure 6: Locations within the cavity where profiles have been extracted 


Prior to quantifying the predictive differences related to the choice of turbulence model, a formal solution verifi- 
cation study was performed to assess the adequacy of the grids used. Both qualitative (simple comparison of results 
on all 3 grids) and quantitative (via the GCI) assessments were performed and are shown in Fig. 7. The quantitative 
assessment based on the GCI is displayed as an error bar attached to the fine grid results, and should be interpreted as a 
bounding estimate of the error due to finite grid resolution. The error assessment was performed for each flow property 
of interest at every other profile station shown in Fig. 6. Only the analysis for the Gatski-EAS model simulations are 
shown in Fig. 7, since this model has the highest degree of non-linearity which typically leads to a larger sensitivity to 
grid resolution. The estimated numerical errors in the streamwise velocity profile are extremely small, suggesting that 
even the coarsest grid considered is sufficient to resolve this flow parameter. The errors extracted for the transverse 
(y-component) velocity predictions are more noticeable due to the small variation of the transverse velocity across the 
domain of interest. The largest errors for this quantity are located in the vicinity of the secondary recirculation zone 
(1st axial station shown). Finally, the only significant source of error associated with the 2"^-order velocity statistics 
is in the shear layer near the front of the cavity. In this region, the turbulence is in a non-equilibrium state as the flow 
transitions from a boundary layer to a free shear layer. 

The solution sensitivity to the turbulence model is shown in the profiles displayed in Fig. 8. The selected stream- 
wise stations are identical to those used for the numerical error assessment. The model-form uncertainty is substantial, 
which was expected given that the region of interest is a large separated fiowfield. Focusing first on the streamwise 
velocity profiles, it is clear that the Gatski-EAS model is predicting a slower shear spreading rate than that predicted 
by the Menter models. The peak reverse flow value (near the lower wall of the cavity) given by the Gatski-EAS model 
is also noticeably smaller. The significant disparity in the transverse velocity at the v = 10.43 inch station is a result 
of the differences in the secondary recirculation zone structure noted in Eig. 5. This flow feature was significantly 
larger in the Gatski-EAS results, which explains the smaller velocities (in magnitude) noted near the cavity floor at 
this station. The transverse velocity values obtained from all 3 models are in closer agreement at the v = 11.14 inch 
station, but the profiles again show appreciable differences further downstream. Relatively small deviations in the 
shear layer spreading rate and/or deflection angle will alter how the shear layer bifurcates at the re-attachment location 
along the aft-wall of the cavity. These details (which are small in velocity magnitude) are more pronounced when 
visualizing the transverse velocity profiles. The bottom two rows of images compare the turbulence kinetic energy 
and dominant shear stress (displayed as a covariance) profiles. In general, the Menter-SST model predicts values for 
these properties that are less than or equal to (in magnitude) those produced by the Menter-BSL model. This trend is 
to be expected since the only difference between these models is the stress limiter present in the SST variant; a feature 
that limits the turbulence levels relative to the BSE model. The Gatski-EAS results show an even greater reduction in 
the turbulence levels within the cavity. It should be noted that algebraic stress models tend to have implied built-in 
stress limiters (relative to their linear counterparts).^^ However, the implied shear stress limiting effect associated with 
the Gatski-EAS model is somewhat less than that explicitly present in the Menter-SST model. Hence, some of the 
reduction in the turbulence predicted by the Gatski-EAS model must be due to the anisotropy in the normal stress 
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Figure 7: Velocity components, turbulence kinetic energy, and x,y velocity covariance profile comparisons on each 
grid (error bars show the GCI on the fine grid results) 
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components predicted by this model. 

The velocity statistics (mean and 2"‘^-order correlations) are compared with the PIV measurements in Fig. 9. The 
mean streamwise velocity comparisons show that both of the Menter models are predicting the proper spreading rate 
of the shear layer. The Gatski-EAS model predictions are also reasonable, but the spreading rate is somewhat slower 
than the measurements indicate. The fact that the spreading rates were predicted correctly (or slightly underpredicted) 
was not expected. Standard RAS models {i.e. models without compressibility corrections) are known to overpredict 
the spreading rates of shear layers with a high convective Mach number. The present flow has a convective Mach 
number of 0.9, which is large enough to expect compressibility effects to be significant. Presumably, the short cavity 
length (relative to the thickness of the approach boundary layer) has limited the impact of this compressibility feature. 
It should also be noted that the simulations predict a core flow velocity that is systematically larger (by approximately 
10%) than the PIV measurements. As pointed out in Ref. 13, this may be a consequence of how the flow was seeded 
for the PIV measurements. The flow was seeded in the boundary layer upstream of the cavity through an angled slot 
injector. The flow rate of the seeded air was intentionally kept small (10 SLPM) to minimize disturbances to the flow. 
Given the thick (^^ 0.25 inches) approach boundary layer, it is quite possible that the outer portion of the boundary 
layer may be insufficiently seeded. The net result would be a measurement that is biased towards lower velocity values. 

Overall, the transverse velocity values compare well with the measurements (particularly the Gatski-EAS results) 
with the most notable exception being the transverse velocity values in the shear layer at some of the downstream axial 
stations. The values within the shear layer are well predicted initially, but start to deviate from the measurements at 
the 10.79 inch station, and this trend progressively worsens at stations further downstream. Note also the difference in 
the predictions within the cavity at the v - 10.43 inch station. As pointed out in the discussion of the flow structure 
from Fig. 5, the Gatski-EAS model predicted a secondary recirculation zone that was notably different than what was 
predicted by the Menter models. The favorable agreement between the Gatski-EAS predictions and the measurements 
at this station suggest that this flow structure is a closer representation of what was measured. 

The velocity variances and covariance comparisons show a strong dependence on the turbulence model. With 
few exceptions, the Gatski-EAS model produced results that were closer to the measured values. All of the models 
predicted the streamwise normal stresses reasonably well, but both of the Menter models consistently overpredicted 
the transverse normal stress component by a significant margin. The Gatski-EAS model allows for Reynolds stress 
anisotropies (a feature not captured by linear eddy viscosity models). Hence, this model is able to predict normal 
stress components that are substantially different from one another. The Reynolds shear stress profiles are also sub- 
stantially better predicted by the Gatski-EAS model. The Menter-SST shear stress values are considerably smaller 
(in magnitude) than the Menter-BSL predictions (a direct consequence of the stress limiter in the SST model), but 
they are still in excess of the measurements by at least a factor of 2. As a final note, the measured velocity variances 
are not asymptoting to zero as the freestream is approached. As pointed out previously in the discussion of the mean 
velocity values, the PIV seed particle count may be insufficient in the outer portion of the shear layer. The fact that 
the measured variance values are not decreasing as the freestream is approached appears to support this hypothesis. If 
this issue is indeed responsible for a deficit in the measured mean velocity values (a deficit as large as 10%), then one 
would expect even larger errors in the higher-order statistics. 

Hybrid RAS/LES Formulation 

The hybrid RAS/LES methodology used in this effort is based on the framework originally developed in Ref. 2, 
with subsequent variants described in Refs. 23 and 24. This framework is designed to enforce a RAS behavior near 
solid surfaces, and switch to an EES behavior in the outer portion of the boundary layer and free shear regions. Hence, 
this formulation can be thought of as a wall-modeled EES approach, where a RAS closure is used as the near-wall 
model. The basic idea is to blend the RAS eddy viscosity value with the EES SGS viscosity, along with any transport 
equation that involves a common RAS and SGS property. In this effort, the Menter-BSL k-co RAS modeP^ was 
blended with the one-equation SGS model of Yoshizawa.^^ The Yoshizawa model involves an evolution equation for 
the SGS turbulence kinetic energy, hence the blended expressions that are appropriate for this model combination are: 

Hybrid RAS/SGS viscosity = {F) [RAS viscosity] -r (1 - E) [SGS viscosity] 

Hybrid RAS/SGS k-equation = {F) [RAS k-equation] -l- (1 - E) [SGS k-equation] (5) 

where E is a blending function that varies between 0 and 1 . Note that the transport equation for the RAS specific 
dissipation rate {co) does not have an SGS counterpart. Hence, the blending is not applied to this equation, and all of 
the terms in this equation that involve the eddy viscosity are evaluated based on the RAS relationships. 

The motivation behind the development of this particular hybrid RAS/LES framework is two-fold. Eirst, the 
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Figure 8: Velocity components, turbulence kinetic energy, and x,y velocity covariance profiles comparisons for each 
turbulence model 
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Figure 9: Comparison of fine grid RAS predictions with measurements: PIV measurements (symbols), Gatski-EAS 
(solid red lines), Menter-BSL (dashed blue lines), Menter-SST (dash-dot green lines) 
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blending of two independent RAS and LES closure models offers the flexibility of having an optimized set of closure 
equations for both RAS and LES modes. The second (and more critical) driving factor was the desire to alleviate 
the difficulties associated with the design of grid topologies that are appropriate for purely grid-dependent blending 
paradigms such as those utilized by the original Detached Eddy Simulation (DES)^^’^^ methodology. The movement 
away from simple grid-dependent blending strategies has become more prevalent in recent years. In fact, even the 
DBS formulation has evolved"^’ to incorporate flow-dependent functions to manage the RAS and LES regimes. 

The particular blending function (F) used in this effort is parameterized by the ratio of the wall distance J to a 
modeled form of the Taylor microscale (T): 


F 

1 



- tank 




^ -1 


A = 


V 


F 

C^poj 


tank ^(0.98) 


( 6 ) 


where k is the von Karman constant (0.41), is 0.09, a is a user-defined model constant, and the factor ra^/z“^(0.98) 
is used to force the balancing position of F (i.e. the position where Krf = to 0.99. The value chosen for a 

provides control over the position where the average LES to RAS transition point (defined sls F = 0.99) occurs. If 
resolved LES content is desired for an attached boundary layer, then this constant should be set such that the transition 
point occurs in the region where the boundary layer wake law starts to deviate from the log law. If the transition point 
is enforced at a lower value that is well within the log law region, a dual log layer appears (an effect sometimes 
referred to as the log layer mismatch). Conversely, if the transition point is enforced at a value that extends 
well into the wake region, then the level of resolved turbulent content will be reduced. Details on a procedure to 
analytically determine the value for a that corresponds to a target value is described in Ref. 23. 


Low-Dissipation Numerical Formulation 

The low-dissipation numerical framework utilized by the cell-centered finite volume VULCAN-CED flow solver 
relies on the blending of a non-dissipative advection operator (2"^, 4^^, or 6^^ order) with a dissipative upwind-biased 
operator. The dissipative advection operators available in the solver include the MUSCL^^ (2"^ or 3^^ order), PPM^^ 
(3''^ order), and WENO^^’^^ (4^^ or 6^^ order) schemes. The blending of each class of operator is controlled via a flow- 
dependent function that attempts to discriminate between inviscid and viscous dominated flow features. A variety of 
sensors have been proposed in the literature for this purpose, two of which are considered here. The Ducros^^ sensor 
distinguishes each flow regime by comparing the magnitude of the vorticity (viscous flow marker) with the velocity 
divergence (inviscid flow marker) via the following functional form: 


Iff = 




(v-y) +(vxy) +e2 


(7) 


where 6, a small number to prevent division by zero in uniform flow regions, is defined as 


6 = 


1 X IQ-^ Voo 

Ay, A^^ 


( 8 ) 


The Larssom^ sensor utilizes the heavy-side function to define a binary (on/off) sensor. 


ilf = 


0, e<\ 

1 , e>\ 


(9) 


where the argument 6 is given by 


e = 




maxifx |V X y| , C 2 [cjWol 


( 10 ) 


This sensor retains the sign on the velocity divergence, so it distinguishes between compression and expansion pro- 
cesses. A positive velocity divergence value (expansion process) always returns a zero value. Hence, a judicious 
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selection of the Ci and C 2 coefficients can lead to a sensor that serves as a shock detector if this is the behavior that is 
desired. An alternative definition of 6 was also considered in this effort that ignores the sign of the velocity divergence, 

i.e. 

(|V • V|) 

^ ^ — r ( 11 ) 

max[Cx |V X y| , C2 [c/ W^J) 

This sensor does not distinguish between expansions and compressions (so it is similar to the Ducros sensor from that 
perspective), while still offering some level of fine-grained control through the selection of the C\ and C 2 coefficients. 
This sensor will be termed the modified Larsson sensor in subsequent sections. Further sensor implementation details 
(smoothing, boundary treatment, etc.) are described in Ref. 8 along with a more detailed description of the dissipative 
and non-dissipative advection schemes. 

Hybrid Reynolds-Averaged / Large Eddy Simulation Results 

The hybrid RAS/LES utilized the low-dissipation numerical framework described previously to encourage the 
development of resolved turbulent content. Both the MUSCL (/c=l/3) and WENO schemes were considered for the 
dissipative advection operator, and the 4^^-order symmetric reconstruction approach was chosen for the non-dissipative 
scheme. The flux limiter used for the MUSCL scheme was the UNO limiter of Suresh and Huynh,^^ which can be 
regarded as a second-order extension of the minmod TVD limiter. The viscous fluxes were evaluated in the same 
manner as that used for the Reynolds-averaged simulations described in the previous section. All of the time-accurate 
hybrid RAS/LES solutions were advanced in time using a dual time-stepping approach that combined a Diagonalized 
Approximate Factorization (DAF) scheme^^ for integration in pseudo-time, with a 3-point backwards finite difference 
approximation for integration in real-time. The values selected for the physical time-step and sub-iteration CFL 
constraint were 0.05 [is and 10.0, respectively. The time-step was chosen based on cell residence time considerations 
to ensure that turbulent structures would traverse less than one grid cell length per time-step {i.e. CFL < 1). The 
sub-iteration process was carried out until the residual error dropped 2 orders of magnitude. This level of convergence 
typically required 7-10 sub-iterations for each physical time-step. The value specified for the parameter a that appears 
in the hybrid RAS/LES blending function (see Eq. 6) was 39.1. This value was determined using the analytical 
procedure described in Ref. 23. 

In order to provide an inflow condition with resolved turbulent content, a recycling/rescaling procedure^^ was 
employed. The particular strategy used here can be described as follows: 

• The velocity, temperature, and density fluctuations are extracted from the desired recycling station (see Figs. 1 
and 3) along with the transported turbulence variables {i.e. k and oj). 

• The extracted flow properties are rescaled based on boundary layer scaling laws. 

• The rescaled fluctuations are added to the mean RAS inflow profile (see Fig. 2), and the rescaled instantaneous 
turbulence variables replace their RAS counterparts. 

• An intermittency function is utilized to prevent recycled fluctuations from corrupting the freestream region. 

• A random span wise shift procedure is used to mitigate the tendency of the large scale structures to persist within 
a fixed spanwise path by the recycling process. 

Further details associated with this procedure are documented elsewhere.^^’^"^ 

The initial state of the hybrid RAS/LES flowfield was defined using the converged RAS results. Moreover, ar- 
tificial fluctuations were added to the approach boundary layer flow to reduce the start-up time required to generate 
resolved turbulent content. The simulations were monitored as a function of time to assess the establishment of a 
statistically stationary state prior to gathering flowfield statistics (or analyzing any instantaneous flow properties). Par- 
ticular emphasis was given to the mass flow error and the integrated friction force time histories. A sample time history 
is shown in Fig. 10. Based on these metrics, a statistically stationary state was deemed to be established after roughly 
25,000 time-steps. This integration period corresponded to 10.6 flow-through times, defined as the time required for 
the freestream particles to traverse the length of the cavity. After this time period, flowfield statistics were gathered for 
at least another 75,000 time-steps (31.8 flow-through times). 

The low dissipation numerical framework utilizes a blending strategy to control where the dissipative and non- 
dissipative operators are invoked. In principle, the non-dissipative operators are desired wherever resolved turbulent 
content is present (boundary layers, free shear layers, etc.) and the dissipative operators are desired in regions where 
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Figure 10: Typical hybrid RAS/LES time history used to assess statistical stationarity 


inviscid flow features are dominant (particularly near shock waves). The three blending sensors described previously 
were evaluated for this purpose. An instantaneous snapshot of each sensor is shown in Fig. 11. The Ducros^^ sensor 
is shown in the top image. This sensor simply compares the magnitude of vorticity with the absolute value of the 
velocity divergence, and is a popular choice used in the research community. This sensor is tagging most of 
the inviscid portions of the flowfleld appropriately. However, it is also tagging a signiflcant fraction of the viscous- 
dominated flow regime (particularly inside the cavity) as a region to be treated by the dissipative operators. The 
random manner in which this sensor is being activated in regions where resolved turbulent structure is present is also 
an undesirable feature. The Larsson^"^ sensor, with the Ci and C 2 coeflicients tuned to detect signiflcant compression 
waves (5.0 and 0.05, respectively), is shown in the middle image of Fig. 11. This blending strategy only invokes 
the dissipative numerical scheme where it is absolutely necessary for stability. However, it is prone to exhibit non- 
physical fluctuations in the inviscid regions of the flow (where there is limited natural viscous effects to provide 
sufficient dissipation). The C\ and C 2 coefficients can be altered to force more of the inviscid region to be tagged for 
the dissipative algorithm, but this sensor is sensitive to the sign of the velocity divergence. Hence, regions of flow 
expansions are inevitably going to be tagged as a non-dissipative flow region. This issue is overcome by desensitizing 
the Larsson sensor to compressions and expansions. The results from this modified form of the Larsson sensor is 
shown in the bottom image of Fig. 11. The behavior of this sensor (with C\ = 5.0 and C 2 = 1.0 x 10“^) is qualitatively 
similar to that of Ducros, but with its undesirable features minimized. Based on these observations, the modified 
Larsson sensor was selected for all of the hybrid RAS/LES results described below. 

A comparison of time-averaged velocity profiles in wall units is given in Fig. 12. The van Driest II transformation 
was used in the post-processing of the computational results shown in this figure to account for compressibility effects. 
Also shown in this image is the hybrid RAS/LES blending function (Eq. 6) to illustrate the location in the boundary 
layer where the simulation transitions (in a time-averaged sense) from RAS to LES. The purpose of this comparison is 
to ensure that the value chosen for a is appropriate, and that the recycling process recovers the proper boundary layer 
scaling behavior. A comparison of the hybrid RAS/LES result with the RAS profile shows evidence of a slight log 
layer mismatch just before the start of the wake region, but overall the profile matches the RAS result reasonably well 
(confirming that a feasible value has been chosen for a). There is some deviation from the theoretical log law result 
that used a value of 0.40 for k with B set to 5.1. However, the deviation appears in both the RAS and hybrid RAS/LES 
results, so it is not a result of the recycling process. 

Instantaneous images of the density gradient magnitude are shown for the coarse, medium, and fine grids in Eig. 13. 
This quantity is similar to a Schlieren image and allows both turbulent structures and inviscid structures (shock waves 
and expansions) to be visualized. As evident in the images shown, even the data extracted from the coarsest grid 
displays a rich array of flow structures, which includes weak eddy shocklets emanating from the movement of large 
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Figure 11: Comparison of sensor formulations controlling the blend of dissipative and non-dissipative schemes: 
Ducros (top), Larsson (middle), modified Larsson (bottom) 


coherent structures in the shear layer and a crisp shear layer re-attachment shock system near the back face of the 
fiameholder. All of the grid levels show evidence of the same qualitative fiow features {i.e. no major feature is missing 
as the grid is coarsened), but the finest grid clearly shows a wider range of scales being resolved. 

Figure 14 shows the structure of the time-averaged fiow inside the cavity fiameholder obtained from the hybrid 
RAS/LES model. The result shown was extracted from the fine grid solution that used the MUSCL scheme with the 
UNO limiter as the dissipative portion of the hybrid numerical scheme. Solutions were also obtained with the 4^^ -order 
WENO scheme used instead of the MUSCL approach. However, the differences in the predictions were found to be 
much smaller than the variability noted with grid resolution (to be discussed later), so these results are omitted for 
brevity. For comparison purposes, the steady-state Gatski-EAS RAS result from Fig. 5 has been repeated here, since 
this model proved to be the best performer of the RAS models considered. The hybrid RAS/LES cavity fiow structure 
shows a secondary recirculation zone that is almost as large as that predicted by the Gatski-EAS model. This was an 
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important flow feature noted in the Gatski-EAS results when comparisons were made with measurements. The size 
and the shape of the primary recirculation zones of both predictions are similar as well. 

Figure 15 shows the variation of several statistics involving the velocity components as the number of samples is 
increased. Each time-step (after a statistically stationary state was obtained) contributed to the sample size, so larger 
sample sizes imply an integration further in time. The profiles shown include the mean streamwise and transverse 
velocity components, turbulence kinetic energy (both the resolved and SGS components), and the resolved and SGS 
portions of the x,y velocity covariance. The profile statistics show very little variation suggesting that 31.8 flow- 
through times (75,000 time-steps) was sufficient to converge the statistics. The largest statistical variability appears 
in the mean transverse velocity predictions, and it is less than 0.4% of the freestream velocity. Note that measures of 
the statistical variability can be misleading if the total integration time is far from being sufficient to have captured 
several periods of the largest physical time scales associated with the flow. An analysis of the autocorrelation of the 
fluctuating velocity field (Ru = r)) can provide information in this regard. The autocorrelation (normalized 

by the local turbulence kinetic energy value) was evaluated at 3 strategic locations: 

• A location in the wake region of the approach boundary layer 

• A location within the shear layer that spans the cavity length 

• A location close to the center of the time-averaged primary recirculation zone 

and the results are shown in the left image of Fig. 16. The data used to generate this image covers 10.6 cavity 
flow-through times (1/3 of the integration time used to gather statistics), and were extracted from the coarse grid 
simulations. By examining the autocorrelations, it is clear that the integration time is sufficient for the boundary and 
shear layer regions of the flow {i.e. the velocity fluctuations become decorrelated after roughly 20 [is). However, 
the integration time shown is not sufficient for the region near the center of the primary recirculation zone, since the 
velocity fluctuations are not yet decorrelated after 625 jds. The largest physically relevant time scale (integral time 
scale) at each location considered can be obtained by integrating the area under the autocorrelation curves. If the 
total integration time is large relative to the integral time scale, then an estimate for the convergence of the statistics can 
be obtained from a finite sample size using the approach outlined in Ref. 38. The autocorrelation data was gathered at 
only a few select points in this effort, so quantified values could not be established for the profile data shown in Fig. 15. 
However, given the general behavior of the points selected, the total integration time used to gather the statistics is 
expected to be sufficient except for the immediate vicinity near the core of the recirculation zones (where the mean 
velocity is extremely small). 

The autocorrelations can also provide information about the range of scales present in the simulation. Two different 
points in the wake region of the approach boundary layer (one near the recycling plane and one just before the cavity 
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Figure 13: Hybrid RAS/LES instantaneous images of density gradient magnitude: coarse grid (top), medium grid 
(middle), fine grid (bottom) 



Figure 14: Streamlines extracted from simulations: hybrid RAS/LES (left), RAS with Gatski-EAS (right) 
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Figure 15: Velocity components, turbulence kinetic energy, and x,y velocity covariance profile comparisons with 
increasing statistical sample size 
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Figure 16: Autocorrelation coefficients at various locations (left), power spectral density at several boundary layer and 
shear layer locations (right) 



Figure 17: Hybrid RAS/LES filter width divided by the Kolmogorov length scale (based on RAS data) 


step) were analyzed, as were two points in the cavity shear layer (before and after the shear layer re-attachment shock 
wave). The right image of Fig. 16 shows the resulting Power Spectral Density (PSD) for the range of frequencies 
that were resolved in the coarse grid simulations. An LES should resolve at least a portion of the inertial range of 
the spectrum, and the slope of this portion of the spectrum is expected to be -5/3. The PSD image shows a distinct 
inertial range captured in the spectrum (covering a frequency range of approximately 1 decade) with the expected -5/3 
slope. Hence, even the coarsest grid considered in this effort appears to have the resolution requirements expected for a 
reasonable LES. Figure 17 displays the ratio of the fine grid filter width to the Kolmogorov dissipation scales (estimated 
from the Gatski-EAS RAS data). This image provides an estimate of the additional grid refinement that would be 
required in order to perform a DNS of this fiowfield. If one ignores the near- wall boundary layer requirements (which 
are treated with a RAS framework), then the grid would require up to 100 times more grid cells in each coordinate 
direction to resolve the smallest turbulence dissipation scales. 

As described previously, the sub-grid scale model errors are considered as part of the numerical uncertainty in this 
effort. If the sub-grid model is performing as intended, then the SGS terms should be modeling the missing physics 
(that can not be resolved by the grid) in a manner that duplicates the statistics of a fully-resolved DNS. The portion 
of the SGS terms that deviate from this intended behavior are taken to be lumped with the truncation error of the 
numerical scheme and evaluated using the GCI metric. This uncertainty is displayed as an uncertainty bar attached to 
the fine grid results in Fig. 18, which shows the infiuence of grid resolution on several statistics involving the velocity 
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components. The estimated grid-related uncertainty in the stream wise velocity predictions are relatively small, with 
the largest uncertainty located early in the shear layer development and within the secondary recirculation zone in 
the cavity. The uncertainty associated with the transverse velocity predictions are considerably larger (relative to the 
mean value), and are most pronounced within the secondary recirculation zone. This flow feature has even longer time 
scales than the primary recirculation zone, and the autocorrelation analysis of a point within the primary recirculation 
zone (see Fig. 16) was found to have an integral time scale 1 to 2 orders of magnitude larger than that of the free shear 
layer. Hence, it is quite possible that a much longer integration time is required to establish highly accurate statistics 
in this portion of the flow (much longer than what can be deemed practical). 

Both the resolved component and the SGS component of the 2"^-order velocity statistics (turbulence kinetic energy 
and velocity covariance) are displayed in the last two images of Fig. 18. The SGS components are the smaller values 
(in absolute value) shown in each image. These plots show that the SGS contribution of each 2"^ -order statistic 
is considerably smaller than the resolved contribution for each grid level except at the first station shown. At this 
station the coarse grid SGS components (and to a lesser extent the medium grid values) are comparable to the resolved 
components. This station is just downstream of the cavity step, and represents a region where the RAS modeled 
boundary layer transitions to a scale-resolved shear layer. This transition process is accelerated as the grid is refined, 
as evidenced by the reduction in the SGS components and the escalation of the resolved portion with added grid 
resolution. As a result, there is a significant variability in the turbulence correlations at the v = 10.43 inch station as 
the grid is refined, resulting in the large uncertainty estimates shown. At stations further downstream, the uncertainty 
is reduced (for the most part) for each of the flow variables. Moreover, at the last two stations the variation between 
the grid levels is reduced and the differentials become more monotonic in sign as the grid is refined. This observation 
provides additional confidence in the uncertainty estimates extracted at these stations. 

The hybrid RAS/LES velocity statistics extracted from the coarse, medium, and fine grid simulations are compared 
with the PIV measurements in Fig. 19. The mean stream wise velocity profiles are compared in the top row of images. 
The simulation results are predicting a considerably thinner shear layer width at the first two axial stations, which is a 
trend also seen in the pure RAS results. The window size used to post-process the raw PIV data was 0.07 x 0.07 inches 
(10 times coarser in each direction than the resolution of the fine grid simulations). This resolution would not be 
capable of resolving a shear layer as thin as predicted in the simulations. Hence, it is possible that the difference noted 
here is simply a matter of implicit Altering introduced by the experimental data processing. With the exception of these 
first two stations, the agreement between the simulation results and measurements is quite good, although the coarse 
grid results tend to underpredict the shear layer width. Larger differences are noted in the mean transverse velocity 
predictions, particularly at the first two axial stations. Recall from Fig. 18 that the uncertainty in the transverse 
velocity predictions was large in this region of the cavity, and that the autocorrelation analysis suggested that this 
may be due to an inadequate integration time. The mean transverse velocities are small and have little variability 
throughout the cavity, so errors from any source will be more noticeable in this parameter. A comparison of transverse 
velocity predictions with the RAS results from Fig. 9 show an improvement over the Menter model predictions, but 
no measurable improvement can be claimed over the Gatski-EAS model results. The 2"^-order velocity statistics are 
compared in the last three rows of images. The simulation values shown include both the resolved and modeled SGS 
components of each correlation. For the most part, the hybrid RAS/LES results are overpredicting the magnitude of 
the 2"^-order velocity statistics. Surprisingly, the Gatski-EAS RAS model predictions almost uniformly outperformed 
the hybrid RAS/LES predictions for these correlations. However, the hybrid RAS/LES results are encouraging from 
the standpoint that the spread in the results obtained from each grid level (which includes the uncertainty associated 
with the performance of the SGS model) is considerably less than the epistemic uncertainty associated with the choice 
of RAS model. Moreover, the error in the hybrid RAS/LES is reducible (given adequate computational resources). 
At least for the latter half of the cavity where the uncertainties were deemed quantiflably trustworthy, the uncertainty 
estimates for the higher order statistics shown in Eig. 18 suggest that finer grids would tend to drive the results in a 
direction closer to the measured values. 
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Figure 18: Velocity components, turbulence kinetic energy, and x,y velocity covariance 
grid (error bars show the GCI on the fine grid results) 
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Figure 19: Comparison of hybrid RAS/LES predictions with measurements: PIV measurements (symbols), fine grid 
(solid red lines), medium grid (dashed blue lines), coarse grid (dash-dot green lines) 
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Summary and Future Work 


Reynolds-averaged and hybrid Reynolds-averaged/large-eddy simulations have been performed for a model scram- 
jet cavity flow experiment. Particle Image Velocimetry data was taken at the centerplane of the cavity, allowing for 
detailed comparisons of the simulated cavity flow structure with measurements. The purpose of the computational 
effort was to assess the state-of-the-art for both RAS and hybrid RAS/LES predictions for this flowfleld, which is a 
relevant one to the scramjet research community. The hybrid RAS/LES assessment was of particular interest, since 
recent development activities associated with the VULCAN-CED code have focused on enhancing these capabilities. 
As a part of the state-of-the-art assessment, the numerical errors were estimated for each simulation approach. 

The Reynolds-averaged simulations of this flow were performed using two linear eddy viscosity models (Menter- 
BSL and Menter-SST k-co models) and one non-linear model (k-oj based explicit algebraic stress model of Gatski and 
Rumsey). The numerical errors were shown to be negligibly small relative to the model-form variability associated 
with the choice of turbulence model. In general, the Menter models did a slightly better job of predicting the mean 
streamwise velocity profiles within the free shear layer that spans the cavity. However, the explicit algebraic stress 
model proved to substantially outperform the linear models when predicting all of the other measured properties 
(mean transverse velocity and velocity variances/covariances). The most prominent feature present in the algebraic 
stress model that is missing in the linear models is the ability to predict normal Reynolds stress anisotropies. These 
features have a pronounced effect in applications that involve secondary motions, so it was no surprise that results 
obtained deviated from the linear eddy viscosity model predictions. However, the drastic improvement noted in the 
predictions for almost every flow feature analyzed was not expected. 

In principle, the hybrid RAS/LES methods have the potential to reduce the uncertainty of RAS (particularly for 
massively separated flows like the cavity flow considered here) by resolving a substantial fraction of the turbulent 
flowfleld. However, to achieve this improvement one has to ensure that the flowfleld is sufficiently grid resolved 
such that the sub-grid scale model (either explicit or implicit) can reasonably be expected to model the unresolved 
scales. Assessing whether this goal has been achieved is no easy task due in part to the difficulty with separating pure 
numerical errors from sub-grid modeling errors. The numerical uncertainty in this effort made no attempt to separate 
the two, and instead both errors were treated as a common error source that is reducible via grid resolution. The 
results of this assessment suggest that the overall simulation uncertainty of the scale-resolving predictions (numerical 
plus sub-grid model) was less than the overall simulation uncertainty (numerical plus turbulence model) of the RAS 
predictions. In general, the hybrid RAS/LES results showed clear improvements over the RAS predictions that utilized 
standard linear eddy viscosity models. However, if one assumes that the measurement uncertainty (which was not 
provided) is small, then the best performing model for this flowfleld (in terms of time-averaged statistics) was the 
explicit algebraic stress RAS model. Euture efforts will consider the effects of fuel injection without the effects of 
combustion. This will introduce an additional modeling issue (the Reynolds mass flux vector) which is known to 
compromise the RAS closure predictions. The introduction of combustion chemistry will also be considered which 
complicates the modeling for both the RAS and hybrid RAS/LES efforts. 
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